!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! Computes various statistics for \c dg-comp.
!!
!! \details
!!
!! Various statistics are computed. The general principle is that this
!! module directly implements the computation of the statistics which
!! are not related to the turbulent model, and uses an external
!! function to compute the turbulence related statistics. By doing
!! this, different turbulent models can specify different statistics.
!!
!! Presently computed statistics are:
!! <ul>
!!  <li> point sections
!!   <ol>
!!    <li> \f$\displaystyle \left<\rho\right>_{\Sigma T} = 
!!     \frac{1}{t-t_0}\int_{t_0}^tdt \frac{1}{|\Sigma|}\int_{\Sigma}
!!     \rho\, d\sigma \f$
!!    <li> \f$\displaystyle \left<\underline{u}\right>_{\Sigma T} = 
!!     \frac{1}{t-t_0}\int_{t_0}^tdt \frac{1}{|\Sigma|}\int_{\Sigma}
!!     \underline{u}\, d\sigma \f$
!!    <li> \f$\displaystyle \left<T\right>_{\Sigma T} = 
!!     \frac{1}{t-t_0}\int_{t_0}^tdt \frac{1}{|\Sigma|}\int_{\Sigma}
!!     T\, d\sigma \f$
!!    <li> \f$\displaystyle \left<\left(\rho-\left<\rho\right>_{\Sigma
!!    T}\right)^2\right>_{\Sigma T} = \left<\rho^2\right>_{\Sigma
!!    T}-\left<\rho\right>_{\Sigma T}^2\f$
!!    <li> \f$\displaystyle
!!    \left<\left(\underline{u}-\left<\underline{u}\right>_{\Sigma
!!    T}\right)\otimes
!!    \left(\underline{u}-\left<\underline{u}\right>_{\Sigma
!!    T}\right)\right>_{\Sigma T} =
!!    \left<\underline{u}\otimes\underline{u}\right>_{\Sigma
!!    T}-\left<\underline{u}\right>_{\Sigma T} \otimes
!!    \left<\underline{u}\right>_{\Sigma T}\f$
!!    <li> \f$\displaystyle \left<\left(T-\left<T\right>_{\Sigma
!!    T}\right)^2\right>_{\Sigma T} = \left<T^2\right>_{\Sigma
!!    T}-\left<T\right>_{\Sigma T}^2\f$
!!   </ol>
!!  <li> side sections: we compute mean values and fluctuations for
!!  all the model equations, and for the energy and momentum equations
!!  we separate viscous and turbulent contributions
!! </ul>
!!
!! \note To compute the time integrals, we define two registers,
!! corresponding to two time levels; see \c update_stats for details.
!<----------------------------------------------------------------------
module mod_dgcomp_statistics

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_sum

 use mod_sympoly, only: &
   mod_sympoly_initialized, &
   ev
 
 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid, t_ddc_grid,  &
   affmap,              &
   el_linear_size

 use mod_sections, only: &
   mod_sections_initialized, &
   t_section_collection,   &
   t_point_section, t_side_section, &
   c_interp_data, &
   validate_section_collection, &
   sync_section_collection

 use mod_atm_refstate, only: &
   mod_atm_refstate_initialized, &
   atm_ref, atm_ref_e, atm_ref_s

 use mod_dgcomp_testcases, only: &
   mod_dgcomp_testcases_initialized, &
   ntrcs, t_phc, phc

 use mod_turb_flux, only: &
   mod_turb_flux_initialized, &
   c_turbmod_progs

 use mod_dgcomp_rhs, only: &
   mod_dgcomp_rhs_initialized, &
   t_bcs_error, &
   dgcomp_tens

 use mod_dgcomp_ode, only: &
   mod_dgcomp_ode_initialized, &
   t_dgcomp_stv, t_dgcomp_ods, t_dgcomp_ode

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dgcomp_statistics_constructor, &
   mod_dgcomp_statistics_destructor,  &
   mod_dgcomp_statistics_initialized, &
   update_stats

 private

!-----------------------------------------------------------------------

 type, extends(c_interp_data) :: t_intp_point
  real(wp), allocatable :: p(:,:)
  real(wp), allocatable :: ref_rho(:)
  real(wp), allocatable :: ref_e(:)
  real(wp), allocatable :: ref_phi(:)
  real(wp), allocatable :: ref_p(:)
  !> work array for the element computations
  real(wp), allocatable :: work(:,:)
  real(wp), allocatable :: pp(:)
  real(wp), allocatable :: tt(:)
 end type t_intp_point
 
 integer, protected :: &
   ir_rm, ir_rf, &
   ir_um, ir_uf, &
   ir_tm, ir_tf

 integer(selected_int_kind(15)) :: nn

 logical, protected ::               &
   mod_dgcomp_statistics_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_dgcomp_statistics'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_dgcomp_statistics_constructor(sc,grid,base)
  type(t_grid), target, intent(in) :: grid
  type(t_base), intent(in) :: base
  type(t_section_collection), intent(inout) :: sc

  integer :: i, ie, np
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized        .eqv..false.) .or. &
       (mod_kinds_initialized           .eqv..false.) .or. &
       (mod_mpi_utils_initialized       .eqv..false.) .or. &
       (mod_sympoly_initialized         .eqv..false.) .or. &
       (mod_base_initialized            .eqv..false.) .or. &
       (mod_grid_initialized            .eqv..false.) .or. &
       (mod_sections_initialized        .eqv..false.) .or. &
       (mod_atm_refstate_initialized    .eqv..false.) .or. &
       (mod_dgcomp_testcases_initialized.eqv..false.) .or. &
       (mod_turb_flux_initialized       .eqv..false.) .or. &
       (mod_dgcomp_rhs_initialized      .eqv..false.) .or. &
       (mod_dgcomp_ode_initialized      .eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dgcomp_statistics_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! define register indexes
   ir_rm = 1
   ir_rf = ir_rm + 1
   ir_um = ir_rf + 1
   ir_uf = ir_um + grid%d
   ir_tm = ir_uf + grid%d**2
   ir_tf = ir_tm + 1

   do i=1,size(sc%sections)
     select type(d=>sc%sections(i)%section_details)
      type is(t_point_section)
       sc%sections(i)%ndata = 2 + grid%d+grid%d**2 + 2
       sc%sections(i)%nregs = 2
      type is(t_side_section)
       sc%sections(i)%ndata = &
         2*(       & ! mean values an fluctuations
          1        & ! density, inviscid
         +2        & ! energy, inviscid and viscous
         +2*grid%d & ! momentum, inviscid and viscous
         +2*ntrcs )  ! tracers, inviscid and viscous
       sc%sections(i)%nregs = 2
     end select
   enddo
   call validate_section_collection(sc)
   do i=1,size(sc%sections)
     sc%sections(i)%regs = 0.0_wp
   enddo
   do ie=1,grid%ne
     if(sc%e2ps(ie)%np.gt.0) then
       allocate(t_intp_point::sc%e2ps(ie)%interp_data)
       select type(d=>sc%e2ps(ie)%interp_data); type is(t_intp_point)
         np = sc%e2ps(ie)%np
         allocate( d%p(base%pk,np),                                &
           d%ref_rho(np), d%ref_e(np), d%ref_phi(np), d%ref_p(np), &
           d%work(2+grid%d,np), d%pp(np), d%tt(np) )
         do i=1,base%pk
           d%p(i,:) = ev(base%p_s(i),sc%e2ps(ie)%xi)
         enddo
         d%ref_rho = matmul(atm_ref(:,ie)%rho,d%p)
         d%ref_e   = matmul(atm_ref(:,ie)%e  ,d%p)
         d%ref_phi = matmul(atm_ref(:,ie)%phi,d%p)
         d%ref_p   = matmul(atm_ref(:,ie)%p  ,d%p)
       end select
     endif
   enddo

   ! Set the time counter
   nn = 0

   mod_dgcomp_statistics_initialized = .true.
 end subroutine mod_dgcomp_statistics_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dgcomp_statistics_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dgcomp_statistics_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_dgcomp_statistics_initialized = .false.
 end subroutine mod_dgcomp_statistics_destructor

!-----------------------------------------------------------------------
 
 subroutine update_stats(sc,t,uuu,ods,ode)
  type(t_dgcomp_stv), intent(in)    :: uuu
  real(wp),           intent(in)    :: t
  type(t_dgcomp_ode), intent(in)    :: ode
  type(t_dgcomp_ods), intent(inout) :: ods
  type(t_section_collection), intent(inout), target :: sc

  integer :: i, ie, id, jd, dd
  real(wp) :: rgcv, theta
  real(wp), pointer :: rnew(:), rold(:)
  real(wp) :: tmp1(0,0,0), tmp2(0)            ! will not be used
  class(c_turbmod_progs), allocatable :: tmp3 ! will not be used
  type(t_bcs_error) :: tmp4                   ! will not beused
  character(len=*), parameter :: &
    this_sub_name = 'update_stats'

   ! Swap the two registers
   do i=1,size(sc%sections)
     sc%sections(i)%regs(:,2) = sc%sections(i)%regs(:,1)
     sc%sections(i)%regs(:,1) = 0.0_wp
   enddo

   rgcv = phc%rgas / phc%cv
   dd = ode%grid%d
   elem_do: do ie=1,ode%grid%ne
     
     if(sc%e2ps(ie)%np.gt.0) then
       select type(d=>sc%e2ps(ie)%interp_data); type is(t_intp_point)
         ! evaluate the solution at the section points
         d%work = matmul( uuu%u(:,:,ie) , d%p )
         d%work(1,:) = d%work(1,:) + d%ref_rho  ! rho
         d%work(2,:) = d%work(2,:) + d%ref_e    ! e
         !d%work(2+1:2+dd,:)                     ! U
         ! Compute the velocity u
         do id=0,dd-1
           d%work(2+1+id,:) = d%work(2+1+id,:)/d%work(1,:)
         enddo
         ! Compute the pressure p
         d%pp = rgcv*(  d%work(2,:)/d%work(1,:)             &
                      - 0.5_wp*sum(d%work(2+1:2+dd,:)**2,1) &
                      - d%ref_phi ) + d%ref_p
         ! Compute the temperature T
         d%tt = d%pp/(phc%rgas*d%work(1,:))

         ! Now the statistics
         do i=1,sc%e2ps(ie)%np
           rnew => sc%e2ps(ie)%s(i)%p%regs(:,1)
           rold => sc%e2ps(ie)%s(i)%p%regs(:,2)
           rnew(ir_rm) = rnew(ir_rm) + d%work(1,i)
           rnew(ir_rf) = rnew(ir_rf) + (d%work(1,i)-rold(ir_rm))**2
           rnew(ir_tm) = rnew(ir_tm) + d%tt(i)
           rnew(ir_tf) = rnew(ir_tf) + (  d%tt(i)  -rold(ir_tm))**2
           do id=0,dd-1
             rnew(ir_um+id) = rnew(ir_um+id) + d%work(2+1+id,i)
             do jd=0,dd-1
               rnew(ir_uf+jd*dd+id) = rnew(ir_uf+jd*dd+id) &
                 + (d%work(2+1+id,i)-rold(ir_um+id)) &
                  *(d%work(2+1+jd,i)-rold(ir_um+jd))
             enddo
           enddo
         enddo
       end select
     endif

   enddo elem_do

   ! Side sections
   call dgcomp_tens( tmp1,tmp2,tmp3,tmp4 , &
          ode%grid, ode%base, ode%bcs, ode%viscous_flow, ode%turbmod, &
                 t, uuu%u, uuu%im, uuu%turbmod_progs, &
                    ods%turbmod_diags, sc=sc )

   ! Now the reduction
   call sync_section_collection(sc,mpi_sum,1,'ps')

   ! Note: the following part could be done by one processor per each
   ! section, with a final broadcast.

   ! Include 1/|Sigma| and combine with previous time steps
   theta = 1.0_wp - real(nn,wp)/real(nn+1,wp) ! time integral
   do i=1,size(sc%sections)
     sc%sections(i)%regs(:,1) =                    & ! space integral
        sc%sections(i)%regs(:,1)/sc%sections(i)%svol
     sc%sections(i)%regs(:,1) = theta *sc%sections(i)%regs(:,1) &
                      + (1.0_wp-theta)*sc%sections(i)%regs(:,2)
   enddo

   ! Update the time step counter
   nn = nn + 1
 
 end subroutine update_stats
 
!-----------------------------------------------------------------------

end module mod_dgcomp_statistics

